/*
NOTES TO DATA ANALYSIS:
	1. Treatment effect regressions by social desirability score (IRT score) for pooled treatments
*/

clear all
set more off
capture log close

********************************************************************************

log using "$logfiles/analysis_itt_sds_`c(current_date)'.log", replace

use "$covidclean/smscovid_clean.dta", clear

***Keep relevant observations: consent taken and above 18 years of age, non-missing treatment arms
	keep if consent==1 & age>=18
	keep if ~missing(treatment_arm)
	count //3964

***HOUSEKEEPING - BEHAVIOR TREATMENTS, OUTCOMES, ROW TITLES
		
	***Create dummies with treatment arms interacted with behavior
		foreach v in $treat_pool $treat_frames $treat_timings $treat_arms {
			gen `v'_sd = `v'*behavior_sd
			gen `v'_hw = `v'*behavior_hw
		}
		
	***Outcomes		
		gen withothers = inlist(withwhom, 3, 8) if !missing(withwhom)
		gen whereother = inlist(whereabout, 2, 3) if !missing(whereabout)
		global othsdhw act_sd contact_3 list_sd contact_c_3 withothers whereother act_hw washhands_1 list_hw washhands_c_1
		global allsd act_sd contact_3 list_sd contact_c_3 withothers whereother
		global allhw act_hw washhands_1 list_hw washhands_c_1
		
	***Row titles
		global rtitle_pooled "Treatment - SD \ "" \ "" \ Treatment - HW \ "" \ "" "
		global rtitle_pooled_sds "Treatment - SD x SDS \ "" \ "" \ Treatment - HW x SDS \ "" \ "" "

		
		
********************************************************************************
***CORRELATION BETWEEN DIFFERENT OUTCOMES

	*** Social distancing outcomes
	
		pwcorr $allsd if treatment_pooled == 0 & interviewdelta == 0, sig
		mat sd_sig = r(sig)
		local sdcount: word count $allsd
		mat sd_pwcorr = r(C)
		*Statistical significance
		mat sd_stars = J(`sdcount',`sdcount',.)
		forval i = 1/`sdcount' {
			forval j = 1/`sdcount' {
				if `i' <= `j' {
					if (sd_sig[`j',`i'] < 0.1)  mat sd_stars[`j',`i'] = 1
					else mat sd_stars[`j',`i'] = 0
					if (sd_sig[`j',`i'] < 0.05) mat sd_stars[`j',`i'] = 2 
					if (sd_sig[`j',`i'] < 0.01) mat sd_stars[`j',`i'] = 3
				}
				mat sd_stars[`i',`i'] = 0	
				if `i' > `j' {
				mat sd_stars[`j',`i'] = 0	
				}
			}
		}
		*Remove all duplicate values in the matrix
		forval i = 1/`sdcount' {
			forval j = 1/`sdcount' {
				if `j' < `i' {
				mat sd_pwcorr[`j',`i'] = .
				}
			}
		}
		*Control means
		foreach v of varlist $allsd {
			mean `v' if treatment_pooled==0 & $sample
			mat mean = r(table)
			local `v'_c_mean : di %9.2f mean[1,1]
		}
		
		*Output 
		global rtitle_sd "Main \ Direct \ List \ Community \ With Others \ At Others"	
		
		global font "basefont(scriptsize) statfont(scriptsize; scriptsize; scriptsize; scriptsize) rtitlfont(scriptsize; scriptsize; scriptsize; scriptsize) ctitlfont(scriptsize) notefont(scriptsize)"
		
		frmttable using "$tables/sd_pwcorr", statmat(sd_pwcorr) store(sd_pwcorr) annotate(sd_stars) asymbol(*,**,***) sdec(3) tex fr replace $font ///
					rtitle($rtitle_sd) ctitle("", "Main", "Direct", "List", "Community", "With Others", "At Others") /// 
					addrows("Control mean", "`act_sd_c_mean'", "`contact_3_c_mean'", "`list_sd_c_mean'", "`contact_c_3_c_mean'", "`withothers_c_mean'", "`whereother_c_mean'") hlines(1{0};1{0};1{0}1;1{0}1)
					
		
	*** Hand washing outcomes
		
		pwcorr $allhw if treatment_pooled == 0 & interviewdelta == 0, sig
		mat hw_sig = r(sig)
		local hwcount: word count $allhw
		mat hw_pwcorr = r(C)
		*Statistical significance
		mat hw_stars = J(`hwcount',`hwcount',.)
		forval i = 1/`hwcount' {
			forval j = 1/`hwcount' {
				if `i' <= `j' {
					if (hw_sig[`j',`i'] < 0.1)  mat hw_stars[`j',`i'] = 1
					else mat hw_stars[`j',`i'] = 0
					if (hw_sig[`j',`i'] < 0.05) mat hw_stars[`j',`i'] = 2 
					if (hw_sig[`j',`i'] < 0.01) mat hw_stars[`j',`i'] = 3
				}
				mat hw_stars[`i',`i'] = 0
				if `i' > `j' {
				mat hw_stars[`j',`i'] = 0
				}
			}
		}
		*Remove all duplicate values in the matrix
		forval i = 1/`hwcount' {
			forval j = 1/`hwcount' {
				if `j' < `i' {
				mat hw_pwcorr[`j',`i'] = .
				}
			}
		}
		*Control means
		foreach v of varlist $allhw {
			mean `v' if treatment_pooled==0 & $sample
			mat mean = r(table)
			local `v'_c_mean : di %9.2f mean[1,1]
		}
		
		*Output 
		global rtitle_hw "Main \ Direct \ List \ Community"		
		
		frmttable using "$tables/hw_pwcorr", statmat(hw_pwcorr) store(hw_pwcorr) annotate(hw_stars) asymbol(*,**,***) sdec(3) tex fr replace $font ///
					rtitle($rtitle_hw) ctitle("", "Main", "Direct", "List", "Community") /// 
					addrows("Control mean", "`act_hw_c_mean'", "`washhands_1_c_mean'", "`list_hw_c_mean'", "`washhands_c_1_c_mean'") hlines(1{0};1{0};1{0}1;1{0}1)
					
********************************************************************************
***MAIN ITT REGRESSIONS
********************************************************************************
	
		qui foreach v of varlist $othsdhw {
			
			********************************************************************
			***BY POOLED TREATMENT
			********************************************************************
			
				*Control mean
				mean `v' if treatment_pooled==0 & $sample
				mat mean = r(table)
				local c_mean : di %9.2f mean[1,1]

				*Regression
				noi reghdfe `v' $treat_pool_sd $treat_pool_hw if $sample & ~missing(sds_irt), a($studycontrols $covariates) vce(robust)
				mat b = r(table)
				local tcount: word count $treat_pool_sd $treat_pool_hw				
				forval i = 1/`tcount' {
					local b_`i' = b[1,`i'] //beta
					local se_`i' = b[2,`i'] //asymptotic s.e.
					local a_p_`i' = b[4,`i'] //asymptotic p-value
					local p_`i' = 0 //exact p-value
				}
				
				local n = e(N)
				local r2_a : di %9.2f e(r2_a)
				
				/*Create matrices to store estimates from Randomization inference 
				mat `v'_p_sds_ri_b = J(`tcount',$nbatch,.)
				mat `v'_p_sds_ri_se = J(`tcount',$nbatch,.)
				mat `v'_p_sds_ri_p = J(`tcount',$nbatch,.)*/
				
				
				*Randomization inference using RI dataset created
				preserve
    				forval batch = 1/$nbatch {
					   use if batch == `batch' using "$rianalysis/ri_data_final", clear
					   gen withothers = inlist(withwhom, 3, 8) if !missing(withwhom)
					   gen whereother = inlist(whereabout, 2, 3) if !missing(whereabout)
    				   noi di "Running ITT-pooled for `v' for batch - `batch' at time: `c(current_time)'"
    				   reghdfe `v' $treat_pool_sd $treat_pool_hw if $sample, a($studycontrols $covariates) vce(robust) 
    				   mat ri_b = r(table)
    				   forval i = 1/`tcount' {
            			   local ri_b_`i' = ri_b[1,`i']
            			   *One-tailed comparison
				           if `ri_b_`i'' >= `b_`i'' local p_`i' = `p_`i'' + 1 
					    }  
				    } 
				restore
				
				*Randomization inference using RI dataset created
				preserve
    				forval batch = 1/$nbatch {
					   use if batch == `batch' using "$rianalysis/ri_data_final", clear
					   *use if batch == `batch' using ri_data_final, clear
					   gen withothers = inlist(withwhom, 3, 8) if !missing(withwhom)
					   gen whereother = inlist(whereabout, 2, 3) if !missing(whereabout)
    				   noi di "Running ITT-pooled for `v' for batch - `batch' at time: `c(current_time)'"
    				   reghdfe `v' $treat_pool_sd $treat_pool_hw if $sample & ~missing(sds_irt), a($studycontrols $covariates) vce(robust) 
    				   mat ri_b = r(table)
    				   forval i = 1/`tcount' {
            			   local ri_b_`i' = ri_b[1,`i']
            			   /*mat `v'_p_sds_ri_b[`i',`batch'] = ri_b[1,`i'] //beta
						   mat `v'_p_sds_ri_se[`i',`batch'] = ri_b[2,`i'] //asymptotic s.e.
						   mat `v'_p_sds_ri_p[`i',`batch'] = ri_b[4,`i'] //asymptotic p-value*/
            			   *Two-tailed comparison
				           if abs(`ri_b_`i'') >= abs(`b_`i'') local p_`i' = `p_`i'' + 1 
					    }  
				    } 
				restore
				
				*Compute Exact p-values 
				forval i = 1/`tcount' {
					local p_`i' = `p_`i''/$nbatch
				}
				
				*Imputing coefficients, asymptotic and exact p-values in the output matrix
				mat `v'_p_sds = J(`tcount',3,.)
				forval i = 1/`tcount' {
					mat `v'_p_sds[`i',1] = round(`b_`i'',0.001)
					mat `v'_p_sds[`i',2] = round(`se_`i'',0.001)
					mat `v'_p_sds[`i',3] = round(`p_`i'',0.001)
				}
				
				*Statistical significance
				mat stars = J(`tcount',3,.)
				forval i = 1/`tcount' {
					mat stars[`i',2] = 0
					mat stars[`i',3] = 0
					if (`a_p_`i'' < 0.1)  mat stars[`i',1] = 1
					else mat stars[`i',1] = 0
					if (`a_p_`i'' < 0.05) mat stars[`i',1] = 2 
					if (`a_p_`i'' < 0.01) mat stars[`i',1] = 3 
				}
				
				if "`v'" == "list_sd" {
					forval i = 1/3 {
						mat `v'_p_sds[2,`i'] = .
						mat stars[2,`i'] = 0
					}
				}
				else if "`v'" == "list_hw" {
					forval i = 1/3 {
						mat `v'_p_sds[1,`i'] = .
						mat stars[1,`i'] = 0
					}
				}
				
				*Save matrix
				matsave `v'_p_sds, replace saving path($tables)
				/*matsave `v'_p_sds_ri_b, replace saving path($tables)
				matsave `v'_p_sds_ri_se, replace saving path($tables)
				matsave `v'_p_sds_ri_p, replace saving path($tables)*/
				
				*Output 
				frmttable, statmat(`v'_p_sds) store(`v'_p_sds) substat(2) annotate(stars) asymbol(*,**,***) rtitle($rtitle_pooled) nodisplay sdec(3,3,3) 
								
		}
		
********************************************************************************
***ITT REGRESSIONS, HETEROGENEITY BY SDS IRT SCORE
********************************************************************************

*Create interaction term for pooled treatment*SDS
foreach var in $treat_pool_sd $treat_pool_hw {
	gen `var'_sds = `var'*sds_irt 
}
global treat_sds treatment_pooled_sd_sds treatment_pooled_hw_sds
	
	
		qui foreach v of varlist $othsdhw {
			
			********************************************************************
			***BY POOLED TREATMENT
			********************************************************************
			
				*Control mean
				mean `v' if treatment_pooled==0 & $sample
				mat mean = r(table)
				local c_mean : di %9.2f mean[1,1]

				*Regression
				reghdfe `v' $treat_pool_sd $treat_pool_hw sds_irt $treat_sds if $sample & ~missing(sds_irt), a($studycontrols $covariates) vce(robust)
		
				mat b = r(table)
				local tcount: word count $treat_sds	
				local j = 4
				forval i = 1/`tcount' {
					local b_`i' = b[1,`j'] //beta
					local se_`i' = b[2,`j'] //asymptotic s.e.
					local a_p_`i' = b[4,`j'] //asymptotic p-value
					local p_`i' = 0 //exact p-value
					local j = `j' + 1
				}
				
				local n : di %9.0fc e(N)
				local r2_a : di %9.2f e(r2_a)
				
				/*Create matrices to store estimates from Randomization inference 
				mat `v'_sds_ri_b = J(`tcount',$nbatch,.)
				mat `v'_sds_ri_se = J(`tcount',$nbatch,.)
				mat `v'_sds_ri_p = J(`tcount',$nbatch,.)*/
				
				*Randomization inference using RI dataset created
				preserve
    				forval batch = 1/$nbatch {
					   use if batch == `batch' using "$rianalysis/ri_data_final", clear
					   *use if batch == `batch' using ri_data_final, clear
					   *Create interaction term for pooled treatment*SDS
					   foreach var in $treat_pool_sd $treat_pool_hw {
							gen `var'_sds = `var'*sds_irt 
					   }
					   gen withothers = inlist(withwhom, 3, 8) if !missing(withwhom)
					   gen whereother = inlist(whereabout, 2, 3) if !missing(whereabout)
    				   noi di "Running hetero analysis by SDS for `v' for batch - `batch' at time: `c(current_time)'"
					   reghdfe `v' $treat_pool_sd $treat_pool_hw sds_irt $treat_sds if $sample & ~missing(sds_irt), a($studycontrols $covariates) vce(robust) 
    				   mat ri_b = r(table)
					   local j = 4
    				   forval i = 1/`tcount' {
            			   local ri_b_`i' = ri_b[1,`j']
            			   /*mat `v'_sds_ri_b[`i',`batch'] = ri_b[1,`i'] //beta
						   mat `v'_sds_ri_se[`i',`batch'] = ri_b[2,`i'] //asymptotic s.e.
						   mat `v'_sds_ri_p[`i',`batch'] = ri_b[4,`i'] //asymptotic p-value*/
            			   *Two-tailed comparison
				           if abs(`ri_b_`i'') >= abs(`b_`i'') local p_`i' = `p_`i'' + 1  
						   local j = `j' + 1
					    }  
				    } 
				restore
				
				*Compute Exact p-values 
				forval i = 1/`tcount' {
					local p_`i' = `p_`i''/$nbatch
				}
				
				*Imputing coefficients, asymptotic and exact p-values in the output matrix
				mat `v'_sds = J(`tcount',3,.)
				forval i = 1/`tcount' {
					mat `v'_sds[`i',1] = round(`b_`i'',0.001)
					mat `v'_sds[`i',2] = round(`se_`i'',0.001)
					mat `v'_sds[`i',3] = round(`p_`i'',0.001)
				}
				
				*Statistical significance
				mat stars = J(`tcount',3,.)
				forval i = 1/`tcount' {
					mat stars[`i',2] = 0
					mat stars[`i',3] = 0
					if (`a_p_`i'' < 0.1)  mat stars[`i',1] = 1
					else mat stars[`i',1] = 0
					if (`a_p_`i'' < 0.05) mat stars[`i',1] = 2 
					if (`a_p_`i'' < 0.01) mat stars[`i',1] = 3 
				}
				
				if "`v'" == "list_sd" {
					forval i = 1/3 {
						mat `v'_sds[2,`i'] = .
						mat stars[2,`i'] = 0
					}
				}
				else if "`v'" == "list_hw" {
					forval i = 1/3 {
						mat `v'_sds[1,`i'] = .
						mat stars[1,`i'] = 0
					}
				}
				
				*Save matrix
				matsave `v'_sds, replace saving path($tables)
				/*matsave `v'_sds_ri_b, replace saving path($tables)
				matsave `v'_sds_ri_se, replace saving path($tables)
				matsave `v'_sds_ri_p, replace saving path($tables)*/
				
				*Output 
				frmttable, statmat(`v'_sds) store(`v'_sds) substat(2) annotate(stars) asymbol(*,**,***) rtitle($rtitle_pooled_sds) nodisplay /// 
				addrows("Adjusted $ R^2$", `r2_a' \ N, "`n'" \ "Control Mean", `c_mean') sdec(3,3,3)	
				
		}


********************************************************************************
***COMPILE OUTPUTS
********************************************************************************

	*First compile results for pooled (w/o interaction), then by interactions
	qui foreach v of varlist $othsdhw {	
		outreg, replay(pooled_oth) merge(`v'_p_sds) store(pooled_oth) nodisplay 
			    
	}
	outreg, replay(pooled_oth) ///
		ct("", "Distancing", "", "", "", "", "", "Handwashing", "", "", "" \ "", "Main", "Direct", "List", "Community", "With Others", "At Others", "Main", "Direct", "List", "Community" \ "\bf{Panel A: Pooled Treatment ITT Regressions}", "", "", "", "", "", "", "", "", "", "") ///
		addrows("", "", "", "", "", "", "", "", "", "", "" \ "\bf{Panel B: Pooled Treatment ITT Regressions by SDS}", "", "", "", "", "", "", "", "", "", "" \ "", "", "", "", "", "", "", "", "", "", "") 
		

	qui foreach v of varlist $othsdhw {	
		outreg, replay(pooled_sds) merge(`v'_sds) store(pooled_sds) nodisplay 
	}	
	
	***Globals for outputs
	global font "basefont(scriptsize) statfont(scriptsize; scriptsize; scriptsize; scriptsize) rtitlfont(scriptsize; scriptsize; scriptsize; scriptsize) ctitlfont(scriptsize) notefont(scriptsize)"
	global lines "hlines(1{0}1;{0};{0}1001;1{0}1)"
		
	***Display all results
	cd "$tables"
	outreg using itt_sds, replay(pooled_oth) append(pooled_sds) replace tex fr multicol(1,2,6;1,8,4;3,1,11;11,1,11) $font note("") $lines
	

log close
exit, clear
